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Thermal fluctuations in non-equilibrium steady states generically lead to power law decay of cor¬ 
relations for conserved quantities. Embedded bodies which constrain fluctuations in turn experience 
fluctuation induced forces. We compute these forces for the simple case of parallel slabs in a driven 
diffusive system. The force falls off with slab separation d as ksT/d (at temperature T, and in all 
spatial dimensions), but can be attractive or repulsive. Unlike the equilibrium Casimir force, the 
force amplitude is non-universal and explicitly depends on dynamics. The techniques introduced can 
be generalized to study pressure and fluctuation induced forces in a broad class of non-equilibrium 
systems. 


External objects immersed in a medium typically mod¬ 
ify the underlying fluctuations and in turn experience 
fluctuation-induced force (FIF) [l|. The textbook ex¬ 
ample is the Casimir force Q arising from quantum 
fluctuations of the electromagnetic field; its thermal ana¬ 
log in critical systems U has also been observed in binary 
liquid mixtures (sj. In both cases, the underlying fluctu¬ 
ations are long-range correlated leading to forces that fall 
off as power laws. In the latter (oil/water mixture) this is 
achieved by tuning to a critical point, while the former is 
a consequence of the massless nature of the photon field. 
Generically in a fluid in equilibrium, correlations (and 
hence FIF) decay exponentially and are insignificant be¬ 
yond a correlation length. 

Non-equilibrium situations provide another route to 
long-range correlated fluctuations: Systems which in 
equilibrium have zero or short-ranged correlations {Ceq 
i5®(x) in s dimensions), quite generically exhibit power 
law correlations {Cneq with conserved dynamics 

when out of equilibrium It is thus natural to inquire 
about the nature (strength and range) of FIF in corre¬ 
sponding non-equilibrium settings (where there is no cor¬ 
responding force in equilibrium). Such forces have indeed 
been explored in a number of circumstances, including 
driven granular fluids , shear flow , and in or¬ 
dinary fluids subject to a temperature gradient [ij, [T^ . 
Here, we explore possibly the simplest (and hence ana¬ 
lytically tractable) example of FIF in a system of diffus¬ 
ing particles which are subject to hard core exclusion, 
commonly referred to as the symmetric simple exclusion 
process (SSEP) [T^. 

The setups examined are: (a) The two dimensional 
system shown in Fig. [T](a); infinite in the y direction and 
connected to two reservoirs at a; = 0 and x = L, with 
densities p{0,y) = pi and p{L,y) = pr, respectively. Two 
slabs, a distance d from each other, span the system along 
the X direction, (b) The three dimensional extension of 
this setup depicted in Fig. [I^b), with the two slabs re¬ 
placed by a tube of square cross section, (c) A generalized 
setup in which the slabs (or tube in three dimensions) of 
length R < L, do not necessarily span the entire system. 


Consider first the two dimensional setup of Fig. [TJa). 
For equal reservoir densities, pi = pr, the system is in 
equilibrium, the pressure is uniform throughout the box 
and there is no average force on the slabs. When the 
reservoir densities are different, the (average) density pro- 
hle varies linearly between the two reservoirs, and there 
is an average diffusive current of particles along the x 
direction. Its magnitude is j = DAp/L, where D the dif¬ 
fusion constant of the particles and Ap = {pi — Pr)- Since 
the average density profile is the same on both sides of 
each slab, naively one would again expect no force be¬ 
tween the two plates. However, we find that the pres- 
ence of non-equilibrium long-range correlations [13l Il4| 
for Pi ^ Pr leads to a force between the two slabs, given 
by (for d<^L) 
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Here ks is the Boltzmann constant, T is the tempera¬ 
ture of the surrounding bath and g{pi,Pr) is a positive 
dimensionless function of order one. Note that the force 
is attractive and when expressed in terms of current, or 
the average density gradient Vp = Ap/L, proportional to 

[Vp) . Here the over line denotes an average over the 
steady-state probability distribution. When the three di¬ 
mensional analogue of the above setup is considered and 
a tube with a square cross section connects the two reser¬ 
voirs (see Fig. [IKb)) the force between two parallel slabs 
has the same form, with the same function g(pi,pr). For 
d L the force still decays as 1/d but with a coefficient 
that is smaller by a factor of 2. 

While the force is attractive for SSEP, it can be repul¬ 
sive in other interacting diffusive systems. Using a sim¬ 
plified model, we argue that this is the case in boundary 
driven antiferromagnetic Ising models with spin conserv¬ 
ing dynamics for a certain regime of parameters. 

Finally, our results suggest that when the slabs or tube 
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is of finite extension, R, the force should behave as 
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Here is a positive function of p;, pr and R, while xq is 
the distance of the slabs from the left reservoir. For hard 
core particles the force is attractive and proportional to 
R^. A similar scaling form, but with opposite sign, is 
expected for the boundary driven antiferromagnetic Ising 
model. 

To derive the above results we use the formalism of 
fluctuating hydrodynamics [l3l - [l5j| . In this approach the 
dynamical equation of motion for the particle density can 
be shown, either through a microscopic derivation (for ex¬ 
ample, see M) or through a phenomenological approach, 
to be 


dtp{xR) + dxJ {x,t) = 0 , (3) 

with a stochastic current 

Jfi {x, t) = -Dd^p (x, t) + ^/a~(p)r]^ (x, t) . (4) 

Here, H is a diffusion coefficient and rjfj^ is an uncorre¬ 
lated white noise vector with components p = 1, • • • , s, 
where s the system dimension. The noise has zero 
mean 77^(x,f) = 0, is uncorrelated rjfj^ {x.,t) r/,y {x.',t') = 
Sfj.,uS {t — t') (5 (x — x'); its variance a{p) = 2DkBTp^n{p) 
satisfying a fluctuation-dissipation condition, where k(p) 
is the compressibility of the gas. For diffusing particles 
subject to hard-core exclusion, D is a constant indepen¬ 
dent of the density p, and a{p) = 2Da^p{l — p) [l^. Il4j. 
Here a is a UV cutoff given by the lattice size and we 
use the standard convention where 0 < p < 1 is dimen¬ 
sionless. For simplicity, in what follows derivations are 
mostly restricted to the two dimensions (Fig. [T] (a)); the 
extension to three dimensions is straightforward and for 
we only quote the final results. 

The density is subject to the boundary conditions 
p(0, y) = Pi and p(L, y) = pr at the reservoris, while the 
normal component of the current must vanish on the two 
slabs. In steady-state the average density density profile 
is given by 'p{x,y) = pi + Apx/L, with j = {DAp/L)x. 

It is important to note that the continuum equations 
are valid in the hydrodynamic limit of a corresponding 
lattice obtained as follows: Consider a (hyper-)cubic sys¬ 
tem of volume L‘^ divided into iV® boxes of size where 
^ is a length scale such that = L. The hydrodynamic 
regime corresponds to first letting ^ oo with L/^ = N, 
and then taking the limit iV —^ oo. Equation [3] is valid 
when the system is rescaled and length is measured in 
units where f —?• 0 and = L (T5 |. 

With this in mind and using ideas similar to Refs. [3, 
[3, [T3, El we write the average pressure to leading order 
in the fluctuations as: 

P{pi^)) = ^ ^"lp(x) ■ (5) 



FIG. 1: The setups studied consist of: (a) A two dimensional 
system, infinite in the y direction is connected to two reservoirs 
at a: = 0 and x = L, with densities p(0, y) = pi and p{L, y) = 
pr, respectively. Two slabs, a distance d from each other, span 
the system along the x direction, (b) The three dimensional 
generalization of the above, with the two slabs replaced by a 
tube of square cross section. 


Here, Sp{x) = p(x) — p(x), and primes henceforth indi¬ 
cate derivatives with respect to the density p. To cal¬ 
culate the force between the plates the pressure has to 
be evaluated on both sides of each slab. The hydrody¬ 
namic procedure described above implies that calcula¬ 
tions have to be carried out using Eq. |3] with the cut¬ 
off and then with length scales rescaled at the end of 
the calculation so that L is finite. In practice this im¬ 
plies that any divergent UV contributions to the pressure 
fluctuations need to be removed from the results of the 
calculation. In particular, in equilibrium and using the 
continuum result Sp{x)Sp{x') = a®p(l — p)S{x — x'), one 
has 5p{xY = a®p(l — p)/^^, and the fluctuations do not 
contribute to the pressure as ^ > oo. 

Clearly, in the setup considered, at any location along 
the wall contributions from P(p(x)) cancel. However, as 
we now show 5p(x)^ varies on the opposing faces of each 
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slab leading to a fluctuation induced force. To evalu¬ 
ate the out of equilibrium fluctuation induced contribu¬ 
tion to pressure, note that for the SSEP the fluctuation- 
dissipation relation, with k{p) = gives 
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To compute 6p{x.)'^ on the faces of the slabs, we evaluate 
the fluctuations along the walls in chambers of size L x d 
and L X oo, respectively. The first corresponds to the 
chamber between the walls, and the second to the semi¬ 
infinite surrounding spaces. 

To evaluate Sp{x, we use standard methods EES, 
expanding the equation of motion to linear order in 5p 
about the steady-state profile. To linear order the current 
is 


Jfi {x, t) = -Ddf,6p (x, t) + a/ct {p{x))r]^ (x, t). (7) 

The dynamical equation is then l inear in Jp so that the 
correlation function C^qx') = Sp(x)dp{x') satisfies a 
Lyapunov equation [17l.ll8l|. After several straightforward 
manipulations this can be brought to the form 

(VxT>Vx + Vx'T>Vx')C'rie9(x,x') 

= -^'5(x-x')V^,cr(p(x')), (8) 

where C„eq(x, x') = C{x, x') — 2^tT(p(x'))i5(x —x') is the 
non-equilibrium part of the correlation function. Using 
the average density profile, 'p{x,y) = pi + Apx/L, the 
above equation reduces to calculating the Green’s func¬ 
tion of a Poisson equation: 

(Vx + V2,)C„e,(x, x') = 25(x - x')(Ap)2aVL2 . (9) 


The boundary conditions are such that Cneq = 0 when 
either x and x' are on the reservoirs (since the density on 
the reservoirs is fixed, 5p = Q identically), while on the 
slabs its normal derivative is zero (no current). To calcu¬ 
late the force, density fluctuations have to be calculated 
on the slabs, e.g. Cneq{x) = Cneq{{x, y = 0}, {x, y = 0}), 
evaluated at the same point x = x' on one of the slabs. 
Using standard Fourier methods one finds 

Cneqix) =^AnSm^ , (10) 


with 



In the limit L, one finds to order L/d 
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Conversely, for d ^ T (indicated by the superscript 1) 
and to leading order in djL 
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The Fourier series with A„ oc (mr) ^ corresponds to a 
parabola. For d L this gives 
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which is in fact the expected behavior of a one¬ 
dimensional SSFP [H, Hi- For d^ L, Eq.[T2] leads to a 
constant contribution, corresponding to the d ^ oo limit, 
and a contribution similar to Cngq(x) with a co-efficient 
that is smaller by 2. Using the hydrodynamic proce¬ 
dure described earlier, we observe that Sp({x, y = 0})^ = 
Cneqix). Namely, only the long-range part of the correla¬ 
tion function contributes to the pressure. 

The fluctuation-induced correction to the pressure in 
Eq.m is the product of two factors: the first (given in 
Eq. ini) is positive, while the second (from Eq. [HI) is neg¬ 
ative. This leads to a negative contribution to pressure, 
corresponding to attraction between the slabs. In the 
limit d <C T the contribution from the semi-infinite sur¬ 
rounding spaces is negligible. Integrating the local pres¬ 
sure over the slab leads to a fluctuation-induced force 


^ = /dxi -P"|p(x)cLg(a:) 

^ kenApf z{l-z) 

d Jo ^2(l-p(z))2’ 


(15) 

(16) 


as proposed in Eq. [T] Here p(z) = pi + Apz. Evaluating 
the integral shows that the total force is a concave func¬ 
tion, vanishing at pi = p^. It is straightforward to use 
the above results to verify that in the limit d L the 
force decays in the same form with a co-efficient that is 
smaller by 2. The calculation can be repeated in three 
dimensions for the configuration depicted in Fig. [T] The 
force is now calculated between two opposite slabs, say in 
the y direction and yields the exact same result as above. 

The negative result in Eq. [H] may appear counterin¬ 
tuitive, since it originates from a computation of 5p{x)^. 
To validate this conclusion, and the underlying hydro- 
dynamic procedure, we performed Monte-Carlo simula¬ 
tions on a two-dimensional square lattice and measured 
the pressure along the slab (see Appendix [A] for details). 
The results in the limit d/L <^1 and for different lattice 
sizes are shown in Fig. [21 The numerics compare well 
with the theoretical predictions. 

Equation [5] suggests that the pressure, and there¬ 
fore the force, can be either positive or negative, de¬ 
pending on the relative signs of P" and Cneq- To ex¬ 
plore this further we carry out a perturbation theory 
in Ap for a general model with a density dependent 
diffusion constant D{p). The equation for the average 
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FIG. 2: Numerical results for the fluctuation induced pressure 
in two-dimensions, given by the integrand of Eg. 1151 multiplied 
by Ld. Here pL = 0.1, Ap = 0.6 and three values of L, d such 
that d L are shown. The units are chosen such that the 
lattice spacing is set to a = 1 and ksT = 1. The solid lines 
depict numerical results, while the dashed line is the analytic 
calculation. The numerical method for measuring the pressure 
is described in Appendix El 


density is then V (D(p(x)) • Vp(x)) = 0, and the Lya¬ 
punov Eq. [5] now has D as a function of p. Setting 
p{x) = Pi + pi{x)Ap -h p 2 {x){ApY + • • •, it is straight¬ 
forward to show that to order (Ap)^ the final result in 
Eq. [14] is replaced by 


the sign of the force. Consider for example a model with 
D = k'^il - q^{p - po)), cr(p) = r2(l -h P{p - po)^) and 
boundary conditions with pi = pQ. While we are not 
aware of a direct microscopic realization of this formula, it 
can be considered as an approximation for an Ising model 
with repulsive interactions evolving under Kawasaki dy¬ 
namics, with p denoting, say, the density of down spins. 
There, it is known that in one dimension a(p) has a mini¬ 
mum around some po which depends on the parameters of 
the model, with D{p) peaked around po [13, (On gen¬ 
eral grounds this behavior is expected to persist in higher 
dimensions.) Using the above expressions it is straightfor¬ 
ward to check that to order (Ap)^ the fluctuation induced 
force, F ~ * , is repulsive. 

The non-extensivity of the force in Eq. [1^] is somewhat 
surprising, and different from say the critical Casimir 
force which behaves as E oc for generalized 

slabs of side L in s dimensions [T|. This is because Cneq 
scales inversely with the volume of the confining box, re¬ 
sulting in a local pressure that vanishes for a large slab. 
As such, we expect this force to be more relevant to small 
inclusions as opposed to macroscopic slabs. While the ex¬ 
act solution of the force between two inclusions is beyond 
the scope of this paper, we can provide an estimate based 
on dimensional grounds. To this end, we consider paral¬ 
lel slabs of dimension R, and neglect the fluctuations of 
density at the open sides of the corresponding enclosure. 
One is then left with evaluating the pressure fluctuations 
in a chamber of size R x with boundary densities 
specified by the mean density at the edges of the slab. It 
is then straightforward to see that in the limit d L the 
force is now given by (for SSEP) 
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resulting in a force 
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irrespective of dimension s, where p{z) = pi + {Ap){zo + 
Rz)/L as advertised in Eq. |21 
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where the derivatives with respect to the density are eval¬ 
uated at pi. The second term on the right-hand-side 
shows the explicit dependence of the results on the dy¬ 
namics through the appearance of the diffusion coeffi¬ 
cient. Moreover, there are no apparent restrictions on 


We thank M. Kolodrubetz and A. Polkovnikov for valu¬ 
able discussions and suggestions. A A and YK are sup¬ 
ported by BSE and ISF grants. MK is supported by the 
NSF through grant No. DMR-12-06323. 


[1] M. Kardar and R. Golestanian, Rev. Mod. Phys. 71, 
1233-1245 (1999). 

[2] H.B.G. Gasimir, Proc. K. Ned. Akad. Wet. 51, 793 
(1948). 

[3] G.L. Klimchitskaya, U. Mohideen and V.M. Mostepa- 
nenko. Rev. Mod. Phys. 81, 1827 (2009). 

[4] M.E. Fisher and P.-G. de Gennes, G. R. Acad. Sci. Ser. 


B 287, 207 (1978). 

[5] A. Gambassi, A. Maciolek, C. Hertlein, U. Nellen, L. 
Helden, C. Bechineer, and S. Dietrich, Phys. Rev. E 80, 
061143 (2009). 

[6] G. Grinstein, D.-H. Lee, and S. Sachdev, Phys. Rev. Lett. 
64, 1927 (1990). 

[7] G. Gattuto, R. Brito, U. Marini Bettolo Marconi, F. Nori, 




































































5 


and R. Soto, Phys. Rev. Lett. 96, 178001 (2006). 

[8] C. Cattuto, R. Brito, U. Marini Bettolo Marconi, F. Nori, 
and R. Soto, Phys. Rev. E 76, 011113 (2007). 

[9] M.R. Shaebani, J. Sarabadani, and D.E. Wolf, Phys. Rev. 
Lett. 108, 198001 (2012). 

[10] H. Wada and S.I. Sasa, Phys. Rev. E 67, 065302(R) 
(2003). 

[11] T. R. Kirkpatrick, J. M. Ortiz de Zarate, and J. V. Sen- 
gers, Phys. Rev. Lett. 110, 235902 (2013). 

[12] T. R. Kirkpatrick, J. M. Ortiz de Zarate, and J. V. Sen- 
gers, Phys. Rev. E 89, 022145 (2014). 

[13] For a recent review see B. Derrida, J. Stat. Mech. P07023 
(2007). 

[14] H. Spohn, J. Phys. A: Math. Gen. 16, 4275 (1983). 

[15] J. M. Ortiz de Zarate and J. V. Sengers, Hydrodynamic 
Fluctuations in Fluids and Fluid Mixtures (Elsevier, Am¬ 
sterdam, 2006). 

[16] J. Tailleur, J. Kurchan, and V. Lecomte, J. Phys. A: 
Math. Theor. 41 505001 (2008). 

[17] C. W. Gardiner, Handbook of stochastic methods for 
physics, chemistry, and the natural sciences. Springer 
(1994) 

[18] A. L. Garcia, M. M. Mansour, G. G. Lie and E. Cementi, 
J. Stat. Phys. 47 209 (1987). 

[19] J. S. Hager, J. Krug, V. Popkov and G. M. Schiitz, Phys. 
Rev. E, 63, 056110 (2001). 

[20] G. Bunin, Y. Kafri and D. Podolsky, J. Stat. Phys. 152, 
112 (2013). 

[21] J. R. Dorfman, T. R. Kirkpatrick and J. V. Sengers, 
Annu. Rev. Phys. Chem. 45, 213 (1994). 

[22] J. Machta, I Oppenheim and 1. Procaccia, Phys. Rev. A 
22, 2809 (1980); T. R. Kirkpatrick, E. G. D. Gohen and 
J. R. Dorfman, Phys. Rev. A 26, 995 (1982); Ronald F. 
Fox, J. Phys. Chem. 86, 2812 (1982); R. Schmitz and 
E.D.G. Cohen, J. Stat. Phys. 39, 285 (1985). 

[23] B. Schmittmann and R. K. P. Zia., Phase Transitions and 
Critical Phenomena vol. 17, ed. C. Domb and J. Lebowitz, 
Academic Press, London (1995). 

[24] B. M. Law, R. W. Gammon, and J. V. Sengers, Phys. 
Rev. Lett. 60 1554 (1988). 

[25] R. Dickman, J. Chem. Phys. 87, 2246 (1987). 


Appendix A: Numerical evaluation of pressure 

To measure pressure in a confined SSEP, we per¬ 
form Monte-Carlo simulation on a square lattice of size 


L X (d -I- 1), with lattice constant set to one. The lattice 
sites are then labelled by (nxjriy) with Ux = 1,2, •• • ,L 
and Uy = 1,2, ■ ■ ■ ,d + 1. The pressure of hard-core par¬ 
ticles is purely entropic, and we measure it by a standard 
method used for evaluating entropic pressure in confined 
polymers [l^: To evaluate pressure at site on a wall, 
we introduce a repulsive potential, V = —ksT log X on 
lattice site -1-1). For the remaining sites on this 

row, (n',^,d -|- 1) with ^ Ux, we set V = oo. The 
pressure is then obtained as 


P(n;) = /'dy ^ (Al) 

JO ^ 

where pd+i{n',j,. A) is the average density on site d + 
1). The integral is performed numerically by discretizing 
A € [0,1] into 20 equally spaced values. Monte-Carlo 
simulations are carried out for each value of and A. 


The Monte-Carlo simulations are carried out using 
standard methods: Each lattice site is either occupied 
or empty. At each Monte-Carlo step a site {nx,ny) is 
chosen at random. If the site is occupied and Ux ^ 1, T, 
Uy ^ d,d+l an attempted move of the particle is made to 
one of its randomly chosen nearest-neighbors (with rate 
1 in arbitrary units), as long as the hard-core constraint 
is not violated. 11 Ux = 0 {ux = L), namely near the 
left (right) reservoir, and the site is empty, a particle is 
added with rate a = -r-^ (<5 = r^). If the site is occu- 

l-pi ^ l-Pr^ 

pied, an attempted removal of the particle is made, with 
equal rate as an attempted move to one of the nearest- 
neighbors. This choice corresponds to setting = 1. 
Finally, if = d or d -I- 1 a move is attempted with rate 
min{l, where AE is the difference in potential 

before and after the move. The results presented in the 
main text (Fig. [2|) were obtained using 8 x 10^^ Monte- 
Carlo sweeps. 




